Multistep Methods for Markovian Event Systems

We consider multistep methods for accelerated trajectory generation in the simulation of Markovian event systems, which is particularly useful in cases where the length of trajectories is large, e.g. when regenerative cycles tend to be long, when we are interested in transient measures over a finite but large time horizon, or when multiple time scales render the system stiff. 1. Markovian Event Systems Markovian models are widespread for modeling stochastic phenomena in a variety of domains. Typically, the models are given in a high-level description such as queueing networks, Petri nets, stochastic automata networks, or sets of coupled chemical reactions, amongst many others. In principle, they can be mapped to the stochastic process level in that they are uniquely defined by an initial probability distribution and a generator matrix. But in practice models tend to be very large. The size of the state space typically increases exponentially with the number of system components or, in other words, the model dimensionality. This effect is known as state space explosion and often causes models to be numerically intractable. One major advantage of simulation is that the state space need not be explicitly enumerated. Thus, a model description that reflects the event system character of the model is well suited, in particular for simulation purposes. In almost all relevant cases the structure of the underlying Markov chain is not arbitrary but state transitions correspond to certain events where similar events essentially have the same effect. Hence, they can be taken as specific discrete event systems [4], which provides a structured model description on an intermediate level of abstraction. For Markovian models the events need not be scheduled and the setting of Markovian event systems is also useful for numerical solution [6]. In order to describe a Markovian event system we have to define its state space and to specify all relevant events that may trigger state transitions. It is necessary to define under which conditions a certain event may occur, how it affects the system state and at which rate it occurs. Diverse formal specifications of Markovian event systems can be found in the literature. Here, we adopt the transition class formalism of [10]. Without loss of generality we assume that the state space is S ⊆ N. All events that trigger state transitions are classified according to their effects which yields transition classes. Formally, a transition University of Bamberg, Germany, E-mail: [email protected] class is a triplet C = (U , u, α) where U ⊆ N is the source state space containing all states in which the event or the corresponding state transition, respectively, is possible, u : U → N is the destination state function giving the new state u(x) ∈ N according to the state transition when the event occurs in state x ∈ U , and α : U → R is the transition rate function giving the rate α(x) ∈ R at which the event or transition occurs in state x ∈ U . Any Markovian model can be uniquely described by a set of such transition classes together with an initial distribution. As a queueing network example consider a d-node tandem network with exponentially distributed service times where arrivals occur only at the first node according to a Poisson process with arrival rate λ. The service rates are denoted by μ1, . . . , μd and the buffer capacities by ν1, . . . , νd. Hence, the different types of transitions are arrivals at node 1, moves from node i to node i+ 1, 0 < i < d and departures from node d. Therefore, d+ 1 transition classes are sufficient: C1 = (U1, u1, α1), where • U1 = {(x1, . . . , xd) ∈ N : x1 < ν1}, • u1 : N → N, x → u1(x) = (x1 + 1, x2, x3, . . . , xd), • α1 : N → R, x → α1(x) = λ; Ci = (Ui, ui, αi), i = 2, . . . , d, where • Ui = {(x1, . . . , xd) ∈ N : xi−1 > 0, xi < νi}, • ui : N → N, x → ui(x) = (x1, . . . , xi−2, xi−1−1, xi+1, xi+1, . . . , xd), • αi : N → R, x → αi(x) = μi−1; Cd+1 = (Ud+1, ud+1, αd+1), where • Ud+1 = {(x1, . . . , xd) ∈ N : xd > 0}, • ud+1 : N → N, x → ud(x) = (x1, . . . , xd−1, xd − 1), • αd+1 : N → R, x → αd(x) = μd; It becomes clear that state-dependent rates can be easily incorporated just by corresponding transition rate functions. Also the state space may be infinite, which is then implicitly given by dropping the restrictions on the source state spaces. Phase-type distributed interarrival and service times can be modeled by properly defined transition classes for any change from one to the next phase. As a chemical reaction set consider the enzyme-catalyzed substrate conversion E + S c1 c2 ES c3 ⇀ E + P (1) where c1, c2, c3 denote associated reaction rate constants such that the corresponding state-dependent reaction rate computes as ci times the number of possible combinations of the required reactants. States of corresponding Markovian models are similarly defined as states of a queueing network, namely by the number of molecules of each species. If we successively number the species E,S,ES, P , a state x = (x1, x2, x3, x4) expresses that there are x1 E-molecules, x2 S-molecules, x3 ES-molecules, and x4 P -molecules. Then the transition classes corresponding to the stoichiometric equation (1) are the following.